Draft version October 20, 2010 

Preprint typeset using L£-T^X style cmulatcapj v. 08/13/06 



THE FORMATION OF LOW-MASS BINARY STAR SYSTEMS VIA TURBULENT FRAGMENTATION 

Stella S. R. Offner 

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA 

AND 

Kaitlin M. Kratter 

Department of Astronomy and Astrophysics, 50 St. George Street, University of Toronto, Toronto, Ontario M5R 3H4, Canada and 
Institute for Theory and Computation, Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA 

AND 

Christopher D. Matzner 

Department of Astronomy and Astrophysics, 50 St. George Street, University of Toronto, Toronto, Ontario M5R 3H4, Canada 

AND 

Mark R. Krumholz 

Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA, 95064, USA 

AND 

Richard I. Klein 

Department of Astronomy, University of California, Berkeley, CA, 94720, USA and 
Lawrence Livermore National Laboratory, Livermore, CA, 94550, USA 
Draft version October 20, 2010 

ABSTRACT 

We characterize the infall rate onto protostellar systems forming in self-gravitating radiation- 
hydrodynamic simulations. Using two dimensionless parameters to determine disks' susceptability 
to gravitational fragmentation, we infer limits on protostellar system multiplicity and the mechanism 
of binary formation. We show that these parameters give robust predictions even in the case of 
marginally resolved protostellar disks. We find that protostellar systems with radiation feedback pre- 
dominately form binaries via turbulent fragmentation, not disk instability, and we predict turbulent 
fragmentation is the dominant channel for binary formation for low-mass stars. We clearly demon- 
strate that systems forming in simulations including radiative feedback have fundamentally different 
parameters than those in purely hydrodynamic simulations. 
Subject headings: stars: formation stars: protostellar disks, stars:low-mass 



1. INTRODUCTION 

Th e study of stell ar multiplicity is a centuries old pur- 
suit (|Mitchell|[r767t ). yet we still lack a comprehensive 
theory explaining the formation of the wide range of ob- 
served binaries and multiple systems. In this paper, we 
investigate the relative importance of two paths for bi- 
nary formation using simulations of a turbulent molec- 
ular cloud destined to form a small cluster of low mass 
stars. 

Although there are numerous pr oposed mechan isms for 
binary formation (see for example [Tphline 2003), we fo- 
cus here on two prominent channels: the fragmentation 
of a turbulent core, and the fragmentation of a gravi- 
tationall y unstable disk. Brief l y, the turbul ent core hy- 
pothesis (|Goodwin et al.l I2004J : fFisrierl I2004D posits that 
turbulent fluctuations within a bound core can produce 
multiple non-linear perturbations in density, which ex- 
ceed the local Jeans mass and collapse faster than the 
background core. Multiple peaks within a given core re- 
sult in a bound binary or multiple stellar system. 

The second h y pothesis, disk fra gmentation 
(|Adams et alJ 119891 : iBonnell fc Bate) I1994D . suggests 
that disks subject to sufficiently strong gravitational 
instability might fragment to form one or more com- 
panions. Although previous studies suggested that 
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this process was limi ted to large mass ratio s y stems , 
rece nt work such as [S tamatell os fc Whitworthl (|2009l ) 
and iKratter et al.l (|2010a[ ) has shown that when disks 
continue to be fed at their outer edges, the companions 
can grow substantially. 

The latter process has sh own promise in its appli- 
cation to high mass stars (|Kratter fc Matzner! 120061 : 
iKrumholz et al.1 l2007t ). where rapid mass accretion 
pushes disks towards instability. One of the goals of this 
paper is to investigate whether the conditions in low mass 
turbulent cores are ever sufficiently violent to induce disk 
fragmentation. 

Observations offer limited guidance as to which mech- 
anism dominates, although several lines of evidence in- 
dicate that turbulent fragmentation must account for at 
least some of the binary population. At wide separa- 
tions bmaries_sho\va very broad eccentricity distribu- 
tion (|Melo et al.ll200ll) and random orientations between 
the bi nary orbita l plane and the spins of the component 
stars (Hale 1994). Similarly, disks around wide T Tauri 
star binaries show significant misalignment between the 
planes of the disks and the binary orbit (IJensen et al.1 
l200i IMonin et al.l l200?l ISchoTz et al.l \2QW \. High ec- 
centricity and orbital misalignment are inconsistent with 
these binaries having been created via the fragmentation 
of a single planar disk (unless they result from three-body 
interactions), but are expected from turbulent fragmen- 
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tation. In contrast, binaries at small separations have 
reduced eccentricities and far less misalignment between 
the orbital plane and the spin or disk orientations of the 
individual components. However, this can plausibly be 
explained via tidal circularization of initially misaligne d 
systems (e.g. iLubow fc Ogilvidl2000t iBate et alJ 12000( 1 . 
so these binaries may have formed by either turbulent or 
disk fragmentation. Moreover, misalignment of the spin- 
orbital plane of systems born in disks may also be due 
to n-body interactions and ejections in a multiple sys- 
tem. Consequently, observations are at present unable 
to answer the question of whether disk fragmentation 
can be responsible for a s ignificant fraction of binaries. 
(See iHowe fe "C larke 2009 for a further discussion of the 
observational evidence on this point.) 

Previous analytic work ((Matzner fc Levin] 120 05). and 
numerical studies of isolated disks have su ggested that 
disks around low-mas s stars will be stable (|Bolev et al.1 
120071 iCai et al] 120081 ) . This work has shown that ap- 
propriate thermal treatment is essential, and in partic- 
ular that heating by the central star dominates at dis- 
tances beyond 10's of AU. Recent numerical simulations 
of low-mass star formation have demonstrated that small 
scale fragmentation leading to high multiplicity systems 
is indeed reduced when radiat ive feedback is included 
(|Offner et al.M2009bt [Bati [2009h 1 . 

Despite these advances, resolving both the large-scale 
turbulence of a molecular cloud forming an ensemble 
of stars and the details of disk structure remains com- 
putationally challenging. In order to resolve the opac- 
ity limit for fragmentation, SPH methods require mass 
resolution of 1O _2 M0, while grid based methods must 
achie ve resolutions down to a few AU (|Goodwin et al.1 
12007( 1. Simulations often circumvent this difficulty by 
modeling indi vidual dense cores with simple initial con- 
ditions (e.g., [ Klcsscn & B urkeriJ l2000t IGoodwin et al.1 
12004 iWalch et al] 12009ft . Such work is predicated on 
knowledge of isolated cores properties and sacrifices the 
interaction of core s with their turbulent environment 
(|Offner et al.ll200"8ft . 

In this paper we examine the simulations of Offner et 
al. (2009b, hereafter OKMK09), which follow the birth 
and evolution of protostars in a turbulent molecular 
cloud, in light of the criterion for protostellar disk frag- 
mentation recently proposed by Kratter et al. (2010, 
hereafter KMKK10). We do this with the goal of draw- 
ing conclusions about aspects of stellar binary formation 
which are not resolved within OKMK09. We proceed 
in two steps: first we test the KMKKI0 stability crite- 
ria against non-isothermal, turbulent simulations using 
the high-resolution runs from OKMK09; then we apply 
these criteria in order to predict the outcome of under- 
resolved disk accretion in the low-resolution OKMK09 
simulations. Ultimately, we can assess the relative im- 
portance of turbulence and disk instability in the creation 
of low mass stellar binaries. 

In the upcoming sections we define our key dimension- 
less parameters, outline our testing scheme, describe our 

1 IBate] H2009T ) includes radiative transfer on numerically resolved 
scales, but not radiation feedback from nuclear or accretion lumi- 
nosity of stars, which Offner et al. (2009) show is roughly an order 
of magnitude larger, fn that respect Bate's simulations represent 
an intermediate case between non-radiative and radiative simula- 
tions. 



numerical methodology, and ultimately draw conclusions 
about the nature of low-mass star formation. 

2. CHARACTERIZING ACCRETION AND DISK 
FRAGMENTATION 

2.1. Dimensionless Parameters 

Classically, disks are thought to fragment and become 
unstable when self-gravity becomes sufficiently strong as 
measured by a low value of Toomre's Q parameter: 

Q = ^0, (i) 

where fid is the disk angular velocity, Ed is the disk sur- 
face density, and c s ,d is the disk sound speed (|Toomrel 
11964( 1. More recent work has shown that disk gas must 
also be able to cool efficiently on ce the process of col- 
lapse begins in order to fragment (|Gammiel [2001). The 
latter criterion is a strong constraint in the inner parts 
of disks where viscous heating dominates but is typically 
satisfied in the outer radii of the irradi ated disks that we 
study here (|Kratter et al.ll2008l [2010Hl . 

Although Q remains a good predictor of disk fragmen- 
tation, it is of more limited use for evaluating disk stabil- 
ity in both large-scale star formation simulations and ob- 
servations because it requires precise knowledge of local 
disk properties, which are difficult to model and measure. 

Young protostellar disks are typically driven unstable 
only when they receive mass fas ter than they can proces s 
it down onto the central star ((Matzner fc Levinl 12005). 
Thus it is useful to parametrize disk fragmentation as a 
function of the infall rate from large scales. KMKK10 
demonstrated that one can predict disk fragmentation 
under idealized conditions using two dimensionless num- 
bers normalized to the infall rate. One, the thermal pa- 
rameter £, compares the disk sound speed to the mass 
accretion rate: 



L s,d 

where M- ln is the infall mas s accre tion rate. This relates 
to the IShakura fc Sunvaevj (|1973( 1 a parameter through 
£ = 3a/ Q. Increasing £ at fixed a tends to cause Q to 
decrease until self- gravity becomes strong (Q ~ 1), at 
which point gravitationally induced spiral modes lead to 
an increase in the effective value of a. Beyond critical 
values of £ and a, however, the disk will fragment; prac- 
tically this occurs when £ ~ 2 and a ~ 2/3. For instance, 
a str i ctly i sothermal simulation fed by a slowly rotating 
IShul ((1977ft solution, in which £ = 0.975, should show 
str ong spiral arms but not fra gment, whereas one fed by 
the lFoster fc Chevalierl (|1993l ) collapse solution, in which 
£ 3> 1 at early times, should fragment. 

A second parameter, T, measures rotation by compar- 
ing the orbital period of the infalling gas to the accretion 
time scale: 

= 4 Mn(j)f n . , 

where is the total mass in the star-disk system, f2k,in 
is the Keplerian angular velocity at the circularization ra- 
dius of the infall, and (j)- m is the average specific angular 
momentum. A large value of T, e.g., 10 , implies that 
the system mass changes significantly in only a few disk 
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orbits, whereas small values, T ~ 10~ 4 — 10~ 3 , describe 
disks with a mass doubling time of hundreds to thousands 
of orbits. The magnitude of T affects the disk's aspect 
ratio (H/Rd ~ (r/^) 1 / 3 ), and therefore the winding of 
spiral arms, and the mass of forming fragments. 

Conducting idealized collapse experiments with 
ORION, a 3-D Adaptive Mesh Refinement (AMR) 
gravito-radiation-hydrodynamicscode, KMKK10 probed 
the fragmentation threshold and final outcome (single, 
binary, or multiple) for isothermal disks in which £ and 
r are constant but numerical resolution improves over 
time. They concluded that any disk with £ > 2 — 3 
will fragment, for values of T typical of low mass star 
formation(10~ 3 T S 10 -2 ), and that increasing T has 
a weak stabilizing effect. 

KMKK10 argued that £ and T capture the essential 
aspects of disk thermodynamics in the context of steady 
accretion, and that other parameters, such as the cool- 
ing time compared to ^jj^, influence fragmentation only 
through their effect on £. This assertion is not tested 
within their isothermal simulations, but in §4 we con- 
firm that it is consistent with the behavior of our highest 
resolution (non-isothermal) runs. We then treat £ and 
r as robust predictors of unresolved disk instability and 
fragmentation in our lower resolution simulations. 

3. METHODS 

3.1. Numerical Methodology 

In this section, we give a brief overview of the simula- 
tions completed in OKMK09 that are the basis for this 
work. The calculations are performed with the ORION 
AMR code. The simulation boxes have initially uniform 
density and periodic boundary conditions. Energy is in- 
jected in the form of small velocity perturbations with 
wave numbers in the range 1 < k < 2 for three cross- 
ing times until a turbulent steady state is achieved. The 
calculation we refer to as RT includes radiative transfer 
in the Flux-Limited Diffusion approximation, while the 
second calculation, NRT, uses a barotropic equation of 
state that has no mechanism for the transport of rada- 
tion. Since radiative cooling is very efficient during the 
initial driving phase, all of the gas remains close to 10 K 
in the RT run, and both calculations begin with a similar 
temperature. 

After three crossing times, self-gravity is turned on and 
this point is considered t — 0. The initial mean density 
is 4.46 x 10~ 20 g cm' 3 , the 3D Mach number is 6.6, and 
the total mass is 185 M Q . Energy continues to be in- 
jected at a constant rate to off set the natural tu rbulent 
decay throughout the run (e.g. IStone et al.lll998j ). Stars 
are inserted as point particles on ce the Jeans criterion i s 
exceeded on the maximum level (jKrumholz et al.l feOCM). 
In the RT calculation, the stars are endowed with a sub- 
grid stellar evolution model, which includes accretion lu- 
minosity down to the stellar surface, Kelvin-Helmholz 
contraction, and nuclear burning (see OKMK09 for a de- 
tailed description). The calculations evolve with gravity 
for one cloud freefall time or 0.315 Myr. 

The RT and NRT runs have a minimum AMR cell size 
of 32 AU such that protostellar disks are only marginally 
resolved with ~ 10 cells. Although the simulations refine 
based upon th e Jeans length, thus preventing artificial 
fragmentation dTruelove et al.|[l997f) , they do not resolve 



the d isk scale height with more than a few cells ([Nelsonl 
120061 ). Nonetheless, OKMK09 demonstrated that the 
protostellar accretion rates are converged to within a fac- 
tor of two by running a resolution study of the first form- 
ing protostar in each with an extra three levels of AMR 
refinement and minimum cells spacing of ~ 4 AU. We 
use the high resolution studies RTC and NRTC to assess 
convergence. 

3.2. Data Analysis 

Our parameterization in terms of £ and T requires that 
we distinguish disk matter from the infall, so that we can 
evaluate the disk sound speed c Si d and the angular mo- 
mentum scale (j) in in addition to the accretion rate M- ln . 
The first step is to define the disk-accretion boundary in 
a robust way within our simulations. Imposing a den- 
sity cutoff is not sufficient as the protostellar cores have 
a range of conditions, which produce diverse disk prop- 
erties. Automatic identification is complicated by disk 
flaring, close companions and turbulent filaments of gas 
feeding the disks. For the RT and NRT calculations we 
adopt a relatively low density threshold of 10~ 16 g cm -3 . 
We estimate the total angular momentum vector of the 
gas, and rotate the coordinate frame so that the net an- 
gular momentum vector is parallel to the z axis. Finally, 
we restrict the vertical disk height in the z direction to 
±5 cells from the disk midplane. The combination of a 
gas density cut with these geometric constraints captures 
the flaring and warping of the disk, while excluding gas 
flowing into the disk. These criteria are not sufficient to 
remove abutting disks of very nearby companions from 
the sample, but these events are few and are easily iden- 
tified from changes in the estimate of the disk radius. In 
the case of a multiple system with a single disk we per- 
form the analysis in the reference frame of the primary. 

Once the disk material is identified, we estimate the 
disk mass, radius, mean temperature, accretion rate and 
angular momentum. The disk radius, r^, is defined by 
the distance between the farthest disk cell and the pro- 
tostellar location. The mean sound speed, c s ,d, is given 
by ^Jk-QT^/ Hp, where /i p = 2.33mn is the mean parti- 
cle mass and is the mean mass- weighted temperature 
averaged over the disk. We find that constructing a vol- 
ume weighted mean changes the resulting temperature 
by at most 15%. The accretion rate, M- ln , is calculated 
by taking the difference between the total star-disk sys- 
tem mass, -M*d, from one time step to the next. In rare 
cases, Mjn may be negative if the disk is perturbed by a 
nearby companion or passes through a shock. 

We find that varying the density threshold by factors 
of two has a 10-15% effect on the disk mass and accretion 
rate but may translate into a 50% difference in T and £, 
since these parameters are implicitly sensitive to the disk 
radius and the amount of included or excluded high an- 
gular momentum material near the disk boundary. The 
RTC and NRTC runs, which have better resolved disk 
structure and sharper disk edges, serve as limits on the 
sensitivity to specific disk properties. We discuss in more 
detail in £ 15.11 how the parameters change when defined 
on larger scales. 

4. RESULTS 
4.1. Final Protostellar System Outcomes 
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Figure Q] shows examples of several systems formed in 
each of the runs (RT, NRT, RTC, and NRTC). The sys- 
tems can be divided into three classes: isolated stars, 
binaries and multiples formed via the fragmentation of a 
turbulent core, and binaries and multiples formed via the 
fragmentation of a disk. Because disk radii are of order 
a few hundred AU, whereas turbulent fragmentation of a 
core can occur on scales up to the Bonnor-Ebert radius 
(0.05 pc = 10 4 AU for n H = 10 5 cm- 3 , T = 10 K), we use 
500 AU as a crude dividing line. 

In Figure[2] we show the separation of all possible pairs 
of stars, bound or unbound, as a function of time in the 
simulations including radiative feedback (RT). All but 
one binary system form outside of 500 AU, but within a 
typical core radius. Upon inspection, the exception re- 
sults not from disk fragmentation but rather from frag- 
mentation of a cold filament or stream feeding a more 
massive protostar. We conclude that fragmentation lead- 
ing to binaries in the RT run is consistent with core, not 
disk, fragmentation. 

Note that the pairs that appear to diverge in Figure [5] 
are unaffiliated stars, not spreading binaries. Pairs with 
separations greater than 0.1 pc are excluded since they 
clearly form in different cores and cannot result from 
turbulent core fragmentation. 

The fragmentation history in the calculations without 
feedback is quite different. Turbulent fragmentation still 
occurs on core scales, but now protostellar disks also 
fragment; see Figure ffl As discu s sed in OKMK09, and 
as predicted bv lMatzner fc Levinl (|2005t l. this additional 
fragmentation would have been suppressed by realistic 
heating of the disk via radiative feedback. In run NRT, 
which lacks such heating, disk fragmentation leads to sev- 
eral large multiple systems and the dynamical ejection of 
a few of the smaller companions. Considering that ra- 
diative feedback is present in low-mass star formation, 
we conclude that turbulent core fragmentation is likely 
responsible for the low mass stellar binary population. 

4.2. Simulation Thermal and Rotational Parameters: £ 

and r 

We find that the thermal parameter, £, as illustrated 
in Figure [3] is either constant or decreasing with time for 
each system. As shown, most protostars begin with mean 
values of ~ 1-3, skirting the regime where fragmentation 
is possible. 

The decline in £ over time primarily reflects that the 
accretion rates decrease as the disk temperatures do 
not increase significantly. This is consistent with ob- 
servations that the mean accretion rate falls as proto- 
stars transistion from the Class to Class II stages 
(|Andre fc Montmerld [l99l lEvans et al.ll2009f) . By the 
time protostars reach the Class II stage, their luminosi- 
ties are dominated by stellar rathe r than accretion lumi- 
nosity dWhite fc Hillenbrandll200l . 

The £ value from the NRTC run is shown in Figure [3] 
for reference. It is apparent that without heating from 
radiation feedback £ reaches a value that is a factor of 10 
and 20 times larger than its RT and RTC counterparts, 
respectively. The RTC £ begins within a factor of 2 of £ 
for the corresponding RT protostar. However, the more 
resolved disk is smaller and has a higher mean tempera- 
ture so that the values differ by a factor of nearly 10 at 
40 kyr. Because £ tx T 3 / 2 , any temperature difference is 



magnified. 

The RT disks cluster around mean temperatures of 30- 
40 K. Since the protostellar luminosity is comprised of 
both a stellar and accretion component, heating contin- 
ues and the disks remain fairly warm even in those cases 
that the accretion rate, M- ln: diminishes to ~ 1O _7 M0 
yr _1 . In contrast, the disk in the NRTC run, which has 
a barotropic EOS, remains close to an average temper- 
ature of 10 K, a factor of ~ 4 smaller than in the cases 
with radiative feedback. 

Figure 0] illustrates that the rotation parameter T, 
drops rapidly in time. For the first ~ 10 years, the 
decline arises from increasing disk masses, as high angu- 
lar momentum material settles onto the disk faster than 
the disk can drain material onto the star. 

At later times, decreasing Y results from declining ac- 
cretion rates as the core gas is depleted. This trend 
can be seen in Figure |4] for a number of protostars that 
have evolved for more than 5 x 10 4 years. The NRT and 
RT r values are more similar and fall within the same 
range of parameter space, because Y depends primarily 
on the turbulent initial conditions on large scales, which 
are mostly unaffected by radiative feedback. Both the Y 
and £ curves exhibit similar shape variation as a function 
of time due to their linear dependence on the accretion 
rate. 

Figure [5] shows T as a function of £ at different times 
for the protostellar disks in each simulation. The values 
fall in a fairly narrow strip of parameter space. Although 
fluctuations in the variables obscure the trend somewhat, 
protostars move from the top right of the plot to the bot- 
tom left. For comparison, we plot the Y — £ tracks for col- 
lapsi ng Bonnor-Ebert spheres initially are overdense by 
10% (jFoster fc Chevalidll993ft . (We assume their disks 
are heated by a factor of 3.5 above the core tempera- 
ture, consistent with the average of the RT runs.) The 
r values for the Bonnor-Ebert spheres are calculated as- 
suming solid-body rotation with rotational parameters 
13 = 0.02 and /3 = 0.08, where /3 = E Iot /E gI£LV . Both 
cases lie within the RT data and mimic the evolution 
towards lower £ and Y with time. 

4.3. Test of Fragmentation Criteria 

With the £ and Y values in hand we can compare the 
system outcomes to those predicted by the idealized mod- 
els of KMKK10. We take advantage of both the NRT and 
RT runs to explore disk outcomes across the fragmenta- 
tion boundary in £ — T space. KMKK10 predict systems 
that fall to the left of the solid line in Figure [5] should be 
stable to disk fragmentation, while systems to the right 
should fragment. This is precisely what we find; disks 
to the right of the line, which are primarily NRT runs, 
undergo disk fragmentation, while those to the left, pri- 
marily RT runs, do not. Note that low resolution often 
prevents the fragments from surviving to form binaries 
due to the sink particle algorithm (see H5.3I) ; a similar 
effect was observed in resolution studies in KMKK10. 
We find that in a few cases, RT disks briefly cross the 
fragmentation line as a result of short duration accretion 
variability. 

Figure |6] shows the final value of Y, £, and multiplicity 
for each protostellar system. Here we define a multiple 
system as one with bound neighbors within 2000 AU. 
The RT cases all fall to the left of the fragmentation 




Fig. 1. — Log gas column density for six different protostellar systems, where the star positions are marked by crosses. The column 
density scale runs from 0.1 gem -2 to 100 gem -2 . The scale of the image and run are indicated. 




Fig. 2. — Pair separation as a function of time in 1 kyr bins for 
pairs in the RT simulation, where each symbol indicates a different 
pair. The dashed line at 500 AU indicates a rough boundary be- 
tween the disk and core scale, where the mean disk size in the RT 
simulation is approximately 500 AU. The large majority of pairs 
have separations above 0.1 pc and are not shown here. The large 
(blue) symbols indicate the first time bin. 



Fig. 3. — The thermal parameter, §, as a function of time for 
the protostellar disks in each simulation. The high-resolution RTC 
(dashed, red) and NRTC (solid, blue) runs are shown in bold. The 
solid line without a symbol corresponds to the same first forming 
star that is depicted by the RTC and NRTC runs. The data is 
averaged over 10 4 yr bins. 

line, and all binaries in this region are formed via core 
fragmentation. In contrast, the NRT systems are nearly 
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Fig. 4. — Rotational parameter, T, as a function of time for the 
protostellar disks in each simulation. The lines, symbols, and data 
bins are the same as in Figure [3] 
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Fig. 5. — T as a function of £ for the protostellar disks in each 
simulation. The solid line is the boundary between fragmenting 
and stable disks found by KMKK10. The dashed line indicates 
where we have extended it to lower T values than explored by 
KMKK10. Bonnor-Ebert (B-E) spheres for = 0.08 (top) and 
P = 0.02 (bottom) are given by the (green) dot-dashed lines. The 
lines, symbols, and data bins are the same as in Figure l3l 

all to the right of the fragmentation line, with the multi- 
ple system having the highest £ values. Although there 
are single NRT systems on both sides of the line, we find 
that those in the fragmenting region have previously or 
are still experiencing fragmentation and particle merg- 
ing. The formation of primarily single systems in the RT 
run is consistent with the obs erved multiplicity fraction 
of low mass stars fLada 2006|) . 

4.4. Toomre 's Q and Disk to Star Mass Ratios 
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Fig. 6. — The final values of T and £ at 1 for the protostellar 
disks in each simulation. The high-resolution RT (red) and NRT 
(blue) runs are shown in bold. The hatched area indicates £ and T 
values prone to disk instability as shown in Figure [5] 

Although we rely primarily on £ and T to predict disk 
instability we also examine both the disk to star mass 
ratio and Toomre's Q parameter a s prox ies for disk in- 
stability. Following iKratter et al.l (|2008[) , we define the 
ratio of the disk mass to total system mass as 



fi 



Md 

M, 



* d 



(4) 



When the disk and star are of comparable mass, or fj, ~ 
0.5, we expect that the di sk will be unstable to gravita- 
tiona l fragmentation (e.g.. lShu et al.lfl990t IKratter et al.l 
2008). Figure shows \i as a function of time. 

To calculate a disk-averaged value for Q we use the 
density weighted temperatures to calculate the disk 
sound speed. We estimate the disk-averaged surface den- 
sity as Sd = Md/(27rr|), where the coefficient 2 assumes 
a 1/r disk column density profile. Because very mas- 
sive disks (jj, ~ 0.5) are especially prone to fragmenta- 
tion, we show the trajectory of disks in Q — [i space in 
Figure [8] for the disks in each of the simulations. Thin 
disks with Q < 1 are unstable to fragmentation. Note 
that thicker disks fragment at lower values of Q (of or- 
der 0.7), while non-axis ymmetric gravitational instabili- 
ties can set in at Q ~ 2 (fGoldreich fc Lvnden-Bell|[l965l : 
iSellwood fc Carlberglll984D . 

At early times, just after the onset of collapse, the pro- 
tostellar mass is small compared to the disk. However, 
this phase is brief, approximately a few 10 3 years, and 
these structures may in fact be underesolved flattened 
envelopes rather than rotationally supported disks. Af- 
ter this early phase, the systems settle into a state where 
the disk is approximately one quarter of the system mass. 

As the core mass is accreted and the protostar grows, \x 
declines. Both T and £ are positively correlated with the 
behavior of fi since a declining disk-system ratio signals 
a declining infall rate. 

As suggested by Figure El disks in the RT and RTC 
simulations have Q > 1 without exception and thus lie in 
the stable part of the Q — fJ, parameter space. In contrast, 
the NRTC disk approaches Q ~ 1 from below, suggesting 
that it is extremely unstable at early times. As expected, 
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t (yr) 

Fig. 7. — The disk to system mass ratio fi as a function of time 
for the protostellar disks in each simulation. The lines, symbols, 
and data bins are the same as in Figure [3] The hatched area 
indicates values of fj, prone to disk instablity. 
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Fig. 8. — Toomre Q parameter as a function of the ratio of disk to 
system mass, fj,, for the protostellar disks in each simulation. The 
high-resolution RTC (dashed) and NRTC (solid) runs are shown in 
bold. The hatched area indicates values of fi and Q prone to disk 
instablity. 

early fragment formation occurs in the NRTC disk during 
the first several kyr of the simulation when Q ~ 0.7. As 
shown in Figure /i generally decreases with time. In 
Figure [3 declining [i corresponds to larger Q values and 
increased disk stability. 

4.5. Semi-Analytic Comparison 

Although we measure disk properties directly in these 
simulations, certain properties remain unresolved at 
lower resolution and are likely affected by artifacts such 
as numerical diffusion. In addition, the previous met- 
rics represent global averages, which may disguise smaller 



scale instability. 

To examine the importance of resolution and global av- 
eraging, we make a second estimate of disk stability using 
simple analytic models to predict radius-dependent disk 
properties, while still relying on the well resolved infall 
rates and the luminosities determined in the simulations. 

Disks are driven unstable when they are fed material 
more rapidly than they can process it at a given tempera- 
ture, i.e. when £ ^ 2, as discussed in §[2] To evaluate £ we 
estimate disk temperatures at characteristic rad i i usin g 
the disk irradiation models of Ma tzner fc Levinl (|2005[ ). 
They find that embedded disks absorb a large fraction of 
the emitted starlight because the infall envelope is opti- 
cally thick to visible wavelengths which are caught and 
re-emitted towards the disk in the infrared. Using a ray 
tracing calculation they find that the flux reaching the 
disk surface is approximately: 

F d = -7—5- where (5) 

./^ = O.lle- 35 , (6) 

and e is the accretion efficiency from the core onto the 
star-disk system. We adopt e = 1/3 here (see discussion 
in i )5.2l) . The temperature of a disk in which irradiation 
dominates over viscous heating, and which is optically 
thick at all relevant wavelengths, will be 

This estimate is realistic, since disks near their fragmen- 
tation threshold are ind eed likely to be optically thick 
(|Matzner fc Levinl l2005. discussion below their eq. 37 ) , 
and bec ause viscous heating i s minimal beyond a few tens 
of AU dKratter et al.ll2010bri . In our evaluation we use 
as calculated in our simulation; although this is some- 
what affected by unresolved dynamics through its depen- 
dence on the accretion rate, we believe the error to be 
small. And, since our goal is to evaluate stability in a way 
that is independent of poorly resolved disk radii, we con- 
sider a definite radius-mass relation, r<j = 200(M >t /MQ) 
AU. This scaling follows from the assumption that the 
disk radius scales with the core radius; if cores share a 
common turbulent Mach number, then their maximum 
disk size is proportional to R corc , although there may 
be large fluctuations around this trend. When cores are 
pressure confined and supported by both thermal pres- 
sure and subsonic turbulence, then core radius, and thus 
disk radius, scales linearly with mass at fixed temper- 
ature. Note that the form of this relation is relatively 
unimportant as disk stability is deter mined by the max- 
imum size to which the disk grows (jMatzner fc Levinl 
2005). We adjust the normalization of this scaling to 
match our high resolution runs. In Fig. [9] we plot £ for 
each star in the RT run as a function of the star's cur- 
rent mass. For comparison we also show £ calculated at a 
constant radius of 50 AU for one of the stars. The overall 
variability of £ is due to variable accretion rates and thus 
stellar luminosities, while the increase with mass is due 
primarily to a decrease in temperature as the character- 
istic radius increases. At low masses, the analytic radii 
are significantly smaller than those in the simulations, 
which are somewhat enlarged due to numerical diffusion; 
in reality disks may be larger and less coherent at early 
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Fig. 9. — The value of §d,ir calculated for the estimated ir- 
radiated disk temperature, y., at a mass dependent radius of 
Td = 200(M*/Mq) AU for each particle (black). For comparison 
we also show the £ trajectory for one of the particles at a fixed 
radius of 50 AU (triangles connected with solid, orange line). The 
hatched region indicates approximate values of £ at which systems 
will become unstable for values of 10~ 3 < T < 5 X 10~ 2 . 

times than analytic models predict. We find that the 
radius normalization must be increased to at least 600 
AU, a very large disk size for these stellar masses, before 
a significant number of £ — M values move into the un- 
stable regime. This result is consistent with the global 
trends of £ and T derived in the previous sections. 

Stabilization out to such large radii is due to strong 
external radiation at early times. Otherwise disks be- 
come unstable outside of a much smaller radius, within 
which cooling times are too long and v iscous heating 
alone suffices to suppress fragmenta tion (|Rafikovl I2005I : 
iMatzner fc L^vnl[2005t lClarkdl2009l) 

5. DISCUSSION 
5.1. Sensitivity to Parameter Definitions 

We use the parameters £ and T as proxies for gravi- 
tational fragmentation, but as we discussed in § 13.21 our 
evaluation of these parameters depends somewhat on our 
ability to discriminate disk from infall. To gauge what 
uncertainty this may cause, we we reevaluate T and £ us- 
ing Td, (j) in , and M- ln averaged on (or within) a sphere 
of radius 1000 AU centered on the protostar, well outside 
the actual disk but close enough that the enclosed mass 
is dominated by disk material. 

Figure [10] shows £ for the first forming protostar in the 
RT and RTC simulations with the parameters estimated 
from the fiducial disk definition and from a sphere with 
1000 AU radius. The £ values for the two different geome- 
tries follow a similar progression and they are generally 
within a factor of two. This confirms that the accretion 
rate on small scales is driven by the larger scale charac- 
teristics of the core and is not particularly sensitive to 
the details of the disk geometry or analysis. The differ- 
ence between the RT and RTC case is preserved using 
the spherical geometry and is determined mainly by the 
different mean temperatures in the two cases (see §3). 

We find similar agreement between the two methods of 
evaluation when this analysis is repeated for T: the RT 
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Fig. 10. — The thermal parameter, £, versus time for the first 
forming protostar in the RT and RTC (bold) calculations using the 
fiducial disk definition (solid) and a 1000 AU sphere (dashed). 

rotational parameter in the disk and sphere geometries 
is within a factor of ^2. The two RTC T values follow a 
similar trend and lie within a factor of 3. 

5.2. Simplifying Assumptions 

The caveats of the numerical methods are discussed in 
detail by OKMK09. We revisit two main caveats here in 
the context of these results. 

Magnetic fields are indisputably important in star for- 
mation. Certainly they play a role in launching out- 
flows, which we discuss below, and also serve as a source 
of pressure that may resist and slow the influence of 
gravity. Methods using 3D ideal magnetohydrodynamics 
(MHD) find that the inclusion of magnetic fields sup- 
presses disk fragmentation in the p arameter space ap- 
plicable to low-mass star for mation (|Price fc Bate! [20071 : 
iHennebelle fc T cvssier 2008) , a result that is complemen- 
tary to the role we find for radiative feedback. These 
simulations suggest that even the formation of a disk 
may be suppressed if the field is strong and the angu- 
lar momentum vector is aligned with the field. However, 
ideal MHD implicitly overestimates the field strength in 
collapsing regions by neglecting reconnection and diffu- 
sion. Methods that include either ambipolar diffusion 
or Ohmic dissipation, but do not include radiative feed- 
back, find that fragmentation may still occur in disks, 
particularly if t he disk rotation rate is sufficiently high 
(jMachida etall 120081: iDuffin fc Pudritzl l2009t ) . In sum, 
magnetic fields most likely reduce disk fragmentation, 
and since we find no disk fragmentation when radia- 
tive feedback is included and both accretion and rotation 
rates are low, our conclusions would be unchanged. 

The greatest uncertainty in our results arises from the 
absence of protostellar outflows, which impact both the 
accretion and luminosity. Observations of starless cores 
indicate that the core mass function shape is similar to 
the stellar IMF but is shifted t o higher masses by a fac- 
tor of three (|Alves et all 120071: lEnoch et all 120081) . This 
difference is generally interpreted as an efficiency fac- 
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tor reflectin g the amount of gas launc hed and entrained 
in outflows (jMatzner fc McKed l2000f ) . although it docs 
not imply a one-to-one correspondence between cores and 
stars. 

By making some simple assumptions we can estimate 
an upper limit on the effect of neglecting outflows. We 
first adopt a constant efficiency factor, e = f/3, and 
assume that the accretion time remains constant (i.e., 
M*/M- m — constant). Under this transformation, a star 
in the calculation with final mass 1M© would instead 
accrete with 0.33M; n and have a final mass of O.33M . 
Since accretion dominates the luminosity for most of the 
simulation (OKMK09), we approximate the total lumi- 
nosity by: 



GM*Mn 



(8) 



where / acc is the fraction of accretion energy that is radi- 
ated away and i?* is the protostellar radius. A power-law 
fit of a one-zone stella r evolution model calibrated to the 
evolutionary tracks of iHosokawa fc Omukail (|2009ft gives 
i?* oc M^M-^ 1 Altogether, the luminosity varies as e 16 . 

From equation 7, the disk temperature Td oc 
(e _0 35 i«/r^) 1 / 4 oc e - 89 , where (.7j n ) is ro ughly indepen- 
dent of e (e.g., iMatzner fc McKed 12000) and the disk 
radius in a given core scales inversely with the mass: 
r d oc (j in ) 2 /Mi,. This gives £ oc e~ - 34 and T oc e. Conse- 
quently, outflow mass loss shifts the disks towards lower 
r and higher £ values. For e = f/3, all but two RT sys- 
tems stay to the left of the fragmentation line indicating 
that most systems remain stable. 

If we relax the assumption that the accretion time is 
constant and instead assume that M- m is constant (as in 
the case of a collapsing isothermal sphere), then £ oc e~ - 9 
while r is independent of e. Again only two previously 
stable systems cross the fragmentation line; although, 
several additional systems lie very close to the line. Re- 
ducing the core efficiency thus has a mildly destabilizing 
affect on the disks. These corrected values represent an 
upper limit on the change in the parameters, so it ap- 
pears unlikely that including outflows in the simulations 
will decrease disk stability significantly and lead to in- 
creased multiplicity. 

We note that additional uncertainty may arise from 
the geometric aspect of outflow cavites, which could af- 
fect the nature of the disk illumination. This caveat is 
implicit in the use of equations (5)-(7). It is also pos- 
sible that the interaction of outflows with nearby fila- 
ments and cores will in fact increase the amount of tur- 
bulent fra gmentation, lo wer the characteristic protostel- 
lar mass (|Li et alj|2010t C. Hansen, 2010, private com- 
munication), and thus reduce the mean luminosity fur- 
ther. Preliminary findings indicate that disk fragmenta- 
tion remains uncommon in calculations where additional 
fragmentation is triggered by outflow interactions with 
turbulent filaments (C. Hansen, 2010, private communi- 
cation) . Instead, the amount of turbulent fragmentation 
is elevated relative to disk fragmentation, which supports 
our thesis that disk fragmentation around low-mass stars 



is rare. 



5.3. Resolution Limitations 



In our calculation we take a conservative approach to 
the representation of fragmentation with sink particles. 
Although N-body gravitational interactions between par- 
ticles may be modeled to sub-grid cell accuracy, the grav- 
itational interaction between the gas and the stars is not 
well modeled when the sep arations are only a few grid 
cells ([Krumholz et al.l 12004). We choose to merge close 
particles rather than follow poorly resolved interactions 
of the particles and gas and thus introduce an undeter- 
mined amount of error into the gas and particle dynam- 
ics. As a result, we forfeit resolution of binary systems 
with separations less than ~ 200 AU in the RT and NRT 
simulations. In the RT case mergers happen relatively 
seldom and the particles acheive at least a brown dwarf 
mass before merging. These are all included in Figure 
[21 where their initital pair separations suggest that the 
formation is real rather than numerical. It is also pos- 
sible that unresolved fragmentation occurs close to the 
primary star, although our RTC study suggests that this 
is unlikely in most cases. 

Nonetheless, low resolution gas fragments may form 
sink particles that would have become thermally sup- 
ported or dispersed at higher resolution. This situa- 
tion is more applicable to the NRT simulations, where 
demonstrably unstable disks are produced when radi- 
ation feedback is absent. Although our sink particle 
method excludes formation in high-velocity gas, for the 
NRT run it is likely that some cold fragmentation could 
be reduced by additional sink criteria for the gravita- 
tional potential and gra vitational boundedness of clumps 
(jFederrath et all l2010j) . The NRTC run does exhibit 
fragmentation at small scales which may be indicative of 
an unresolved binary. Thus, our results are qualitatively 
supported by both higher resolution runs and analysis of 
r - £ and Q. 

5.4. Comparison with Observations and Previous Work 

These results affirm previous analytic and numerical 
work modeling protostellar accretion disks. Our simula- 
tion sample is comprised of low-mass stars with accretion 
rates below that required for disk fragmentation as mea- 
sured through £ and T. The most massive RT star is 
~ 2.5 Mq, which is marginally in the reg ime where in- 
stabili ty might be expected according to iKratter et al.1 
(2008). However, they find that such a system would 
likely not become unstable until 0.1-0.2 Myr or later. 
Lower mass stars are expected to fragment at even later 
times or not at all, consistent with the small amount of 
time that the RT stars spend to the right of the fragmen- 
tation line. In contrast, the NRT and NRTC cases are 
highly discrepant with the model predictions. It is clear 
that radiative heating from the primary stars is a crucial 
component in simulations. 

Most of the core fragmentation occurs around separa- 
tions of ~ 0.01 — 0.02 pc. Observations of dense starless 
cores in Perseus suggest that the initital density pro- 
file is relati vely smooth outside the beam resolution of 
~ 1200 AU (jSchnee et al.ll2010h . Two of the If cores in 
the sample are found to have elongation on scales of a 
few thousand AU that may be indicative of unresolved 
fragmentation. A very smooth density distribution could 
preclude wide fragmentation as a means for forming bi- 
nary systems. However, our turbulence-produced frag- 
ments are rare and would be missed at the resolution of 
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iSchnee et~a l. (2010). It is also possible that the observed 
starless cores may never form protostars or are younger 
than the objects we focus on in this paper and may un- 
d ergo fragmentation s ometime in the future. 

iMaurv et al.l (|201Q[ ) probe the multiplicity of Class 
objects with sensitivity down to 50 AU separations. Of 
the five sources in their sample, only one shows evidence 
of a companion. This potential binary has a separa- 
tion of ~ 1900 AU, which suggests that it did not result 
from disk fra g menta tion. Their findings, together with 
lLoonev et al.l (|2000[) . are consistent with our prediction 
that low-mass protostars are not members of high-order 
multiple systems and that any companions mostly likely 
originate from turbulent core fragmentation and form 
with a large initital separation. 

Low-mass star systems like those modeled here are ob- 
served to ha ve a much lo wer binary fraction than higher 
mass stars |Lada 2006j). In particular, only 30% of 
M dwarfs, the most common stellar type, have lower 
mass companions. However, both simulations and ana- 
lytic work predict that the more massive disks around 
high-mass stars undergo frag mentation leading t o ad- 
ditional stellar comp anions (jKratter fc Matznerl 120061 : 
iKrumholz et al.ll2007l ). A significant fraction of the dif- 
ference may be due simply to the absence of disk frag- 
mentation around low mass stars. OKMK09 report a 
single star fraction of 0.8 + 0.2/ — 0.4 or 0.5 ± 3 assuming 
all mergers result instead in a binary. These are likely 
upper and lower limits, respectively, on the actual single 
star fraction, since the distributio n may evolve over tim e 
due to interactions between stars (|Duchene et al.ll2007l ). 
Although the statistics are poor, these v alues are co nsis- 
tent with the sing le star fraction of 0.7 (|Ladall2006[ ). 

Using the simulations we can make predictions for 
the initial range of binary separations. Excluding bi- 
naries below the simulation resolution, these are in fact 
much larger than the typical separation of ~ 50 AU, al- 
though there is a wide distribution around the mean, 
which varies with stellar mass and possibly from re - 
gion to region (jDuauennov fc Mavorlll99ll: lFisher|[200l ; 
? find a gaussian distribution for periods peaked at 
10 5 03 days, implying 48 AU separations for a total 
system mass of 1.5M0. This suggests that significant 
evolution of the separations must take place after the 
core fragments. The discrepancy is even larger for 
very low-mass binaries, whose typical maximum sepa- 
rations are 145q(M to t/M ( T 1 ) 2 AU for systems < O.6M 
(jBurgasser et all 120031; iLafreniere et alll2008D. alt hough 
wide binary exceptions do exist (jRadigan et all I2008L 
2009). We favor a scenario in which close low-mass bina- 
ries form via interaction with the ongoing infall and be- 
tween the two accretion disks; loosely bound companions 
can be stripped by close encounters within the protoclus- 



ter. (The initial subvirial velocity dispersion of the stars 
predisposes young clusters to dynamical interactions that 
will impact the fina l binary distribution; iBate et aTll2003l ; 
lOffner et al1l2009ai 1 

6. CONCLUSIONS 

In this work, we revisit the turbulent, self- 
gravitating radiation-hydrodynamic calculations per- 
formed by Offner et al. (2009) in order to characterize 
the protostellar systems in terms of the thermal param- 
eter, £, and the rotational parameter, T. We first use 
the high-resolution simulations to confirm the stability 
criteria derived in KMKK10 under non-turbulent, ideal- 
ized conditions. The high-resolution simulations further 
demonstrate that the dimensionless parameters are con- 
verged at lower resolution and, thus, they are sufficient to 
describe the evolutionary state of the protostellar disks 
even in the case where the scale height and disk struc- 
ture are not well resolved and where the infall rate is 
fluctuating and turbulent. 

As expected, we find that £ and T are distinct in the 
cases with and without radiative feedback. In the for- 
mer case, fragmentation occurs preferentially on scales of 
^1000 AU rather than within disks. These two param- 
eters indicate that protostellar accretion disks around 
low-mass protostars will be stable against gravitational 
fragmentation for nearly all times. As an independent 
test, we perform a semi-analytic analysis using the sim- 
ulated accretion rates and luminosities, which confirms 
that these disks are stable out to large radii. 

Under a scenario of stable accretion disks, low mass bi- 
nary systems arise as a result of multiple collapse events 
in a turbulent core. It follows that the observed mul- 
tiplicity of low-mass stellar systems arises predominatly 
due to turbulent fragmentation in the parent core. 
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